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ABSTRACT 

We consider the emission of photons from the inner parts of a relativistically expanding plasma out- 
flow, characterized by a constant Lorentz factor, V . Photons that are injected in regions of high optical 
depth are advected with the flow until they escape at the photosphere. Due to multiple scattering 
below the photosphere, the locally emerging comoving photon distribution is thermal. However, as 
an observer sees simultaneously photons emitted from different angles, hence with different Doppler 
boosting, the observed spectrum is a multi-color black-body. We calculate here the properties of the 
observed spectrum at different observed times. Due to the strong dependence of the photospheric ra- 
dius on the angle to the line of sight, for parameters characterizing gamma-ray bursts (GRBs) thermal 
photons are seen up to tens of seconds following the termination of the inner engine. At late times, 
following the inner engine termination, both the number flux and energy flux of the thermal spectrum 
decay as F oc t~ 2 . At these times, the multicolor black body emission results in a power law at low 
energies (below the thermal peak), with power law index F v cx v 3 . This result is remarkably similar 
to the average value of the low energy spectral slope index ( "a" ) seen in fitting the spectra of large 
GRB sample. 



Subject headings: gamma rays:theory — plasmas- 
scattering — X-rays : bursts 

1. INTRODUCTION 

Relativistic outflows, often in the forms of jets, are a 
common phenom ena in many astronomical objects, such 
as microquasars ( |Mirabel fc Rodriguez |[l~994| i Hjellming 
fc Rupenll 1995|), active galactic nuclei (AUNs; |Lmd &; 



B landlord |1985[|Go pal-Krishna et al. 2006 ) and gam ma- 
ray bursts (GRBs; iPaczynski 1986; Goodman 1986). In 



many of these objects, the density at the base of the flow 
is sufficiently high, so that the optical depth to Thomson 
scattering by the baryon-related electrons is much larger 
than unity. As a result, if a source of photons exists deep 
enough in the flow, the emerging spectrum is inevitably 
thermal or quasi-thermal (a Wien spectrum could also 
emerge if the number of photons is conserved by the ra- 
diative processes). These photons escape the flow once 
they decouple fro m the plasma, at the photosphere (e.g., 
Paczyriski J l990|>. 



-radiation mechanisms:thermal — radiative transfer- 



tegrated spectrum is non-thermal. 

Analysis of this effect begins with the non-trivial shape 
of the photosphere in a relativistically expanding plasma 
source. By definition, the photosphere is a surface in 
space which fulfills the following requirement: the optical 
depth to scattering a photon originating from a point on 
this surface and reaching the observer is equal to unity. 
Although it is mathematically a two-dimensional surface 
in space, for spherically symmetric wind, as is considered 
here, this surface is symmetric with respect to rotation 
around the axis to the line of sight. It is therefore appro- 
priate to refer to it as the "photospheric radius" , which is 
a function of the angle to the line of sight, r p h = r p h(0). 

Calculation of the photospheric radius for the scenario 
of a steady, spherically symmetric, relativistic wind char- 
ac terized by a constan t Lore ntz factor T was first c arried 

and later extended by Pe'er 



As a result of the relativistic expansion of the source, 
even if the emitted spectrum (in the plasma comoving 
frame) is purely thermal - i.e., it is not modified by any 
additional, non-thermal radiative process - still the ob- 
served spectrum is not expected to follow the Planck dis- 
tribution. This results from the aberration of light in the 
expanding plasma, and is of pure geometrical nature. It 
is therefore an inherent property of any relativistically 
expanding photon emitting source. The origin of this ef- 
fect lies in the fact that at any given instance, an observer 
sees simultaneously photons that emerge the expanding 
plasma from a range of radii and angles. Therefore, each 
(thermal) photon has its own comoving energy, and is 
seen with a particular Doppler shift. The resulting, in- 
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(2008). For relativistic winds, r 3> 1, and small angles to 
the line of sight, 9 <C 1, the calculation results in a sim- 
ple, yet non-trivial dependence of the photospheric radius 
on the angle to the line of sight: r p h{9) cx (r~ 2 + 2 /3), 
where the proportionality constant depends on the prop- 
erties of the outflow. It is thus found that in relativistic 
expanding wind, the photospheric radius has a strong de- 
pendence on the angle to the line of sight, a fact which 
leads to severa l non-tr i vial c onsequences: for example, it 
was shown by Pe'er | (2008), that for parameters char- 
acterizing GRBs, photospheric emission can be expected 
up to tens of seconds (albeit with a decreasing flux; see 
further discussion below). 

By definition, the photospheric radius provides only a 
first order approximation to the last scattering position 
(=decoupling position) of the thermal photons. This is 
because, in principle, photons can be scattered at any 
point in space in which electrons exist. Therefore, a full 
description of the last scattering position and scattering 
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angle can only be done in terms of probability density 
function P(r, 0). As was shown by Pe'er (20081 and will 
be further discussed here, use of this function is essential 
in calculating the spectrum and flux. The shape of the 
observed spectrum, for example, depends both on the 
individual Doppler shifts of the observed photons, but 
also on the number of photons that undergo the particu- 
lar Doppler shift, an information that is only held in the 
probability density function. 

The probability density function can be considered as 
an extension of the standard use of the photospheric ra- 
dius. Instead of considering a surface in space from which 
thermal photons emerge, one considers the entire space, 
weighted by the finite probability of a photon to emerge 
from an arbitrar y radius r and arbitrary angle 9. Us- 

(2008) calculated the expected 
lie observed thermal photon flux 



ing this function, Pe'er 
temporal decay laws ot t 



and average temperature, following an abrupt termina- 
tion of the inner engine. It was shown there, that at late 
times the thermal flux is expected to decay as F oc t~ 2 
and the average temperature decays as < T >oc t~ 2 / 3 . 
The agreement found between these theoretical predic- 
tions and the late time decay laws of the peak energy 



and flux observed during GRB p rompt emission (Ryde 



2004 2005 Ryde & Pe'er 2009), is one of the key mo 



ti vat ions m studying the properties of thermal emission 
in the context of the GRB prompt emission phase. 

On a more general ground, thermal emission may be 
crucial in understanding the nature of GRB prompt emis- 
sion. In spite of being studied for nearly two decades 
now, the origin of the prompt emission in GRBs is still 
puzzling. In recent years it became clear t hat syn 



chrotron emission, the leading emission mode l ( Rees fc 
Meszaros ||1994[ |Tavani || 19961 ICohen et al. || 19H7| jfj, " 



an 



& Piran 1997; Panaitescu et al. II1999D cannot account 
, II r I ,. II r , , 



for the steepness ot the low energy spectral slopes seen 


dCrider et al. |1997[ 


Preece et al. 1998, 2002 IGhirlanda 


et al. 1 2003| (see, n 


owever 


Bosnjak et al. 2009). 'i'his 



through h eated cloud dPermer fc Bottcher ||2000[ ), jitter 
radiation (j Mcdvcdcv 20001} or decaying magnetic field 
Moreover, in order to account for 



( |Pe'er fc Zhang |2 006 ) . 
the high efficiency ot ll 



high efficien cy ot the prompt e mission seen ( Zhang et 
al. |2007|[Nysewander et al. |2009[ ) using the synchrotron 



model, a highly efficient energy dissipation is required 
which is difficult to be accounted for in the classical in- 
ternal shocks scenario (|Kobayashi et al. |1997 Daigne fe 



MochkovitcTri|1998[ |Lazzati et al ||1999l |Guetta et al. 
2001 p . Contribution from thermal emission thus seems a 



natural way of overcoming both these issues. First , pho- 
tospheric emission is inherent to the fireball model ( Eich- 
ler & Levinson 2000 ; Meszaros & Rees 2000 ; Meszaros et 



al. |2002j IDaigne fc Mochkovitch |2002[|Rees fc Meszaros] 
2005[ ). Second, as it does not originate from any internal 



dissipation, contribution f rom thermal photon s reduces 
the efficiency requirement ( |Ryde fc Pe'er | ;2009). Finally, 
as will be shown here, photospheric emission is able to 
produce low energy spectral slopes which are consistent 
with those observed. 

In spite of the success of the thermal emission model in 
reproducing the late time decay of the peak e nergy and 
flux seen in many GRB's ( Ryde fc Pe'er |2009[ ), the sim- 
ple version of the model suffers several drawbacks. One 



issue that is often raised, is that the low energy spectral 
slopes seen during the prompt phase of many GRBs are 
too shallow to be accounted for by the Rayle igh- Jeans 
tail of the thermal spectrum (e.g., Bellm |2010 ). Indeed, 
the Rayleigh- Jeans tail implies a spectral slope F v oc v 2 , 
while GRB observations show that on the average, the 
low energy photon index in the "Band" function fits is 



al. 



_ -1 ( i .e., F y cx v°- 
|[2006l [2008| . 



sec 



Preece et al. |200"0 Kaneko et 



As discussed above, the observed spectrum resulting 
from a photospheric emission in relativistically expand- 
ing plasma does not necessarily need to be a pure black- 
body, but should, in general, be modified. Therefore, 
a pure blackbody spectrum will in many cases not be 
able to fit the observed spectrum. The main purpose of 
this paper is to calculate the observed spectrum result- 
ing from photospheric emission in a scenario of a steady, 
relativistic outflow. In fact, as we show below, at late 
times following the termination of the inner engine, the 
resulting spectrum is expected to be close to a power law 
below the thermal peak, with power law index F v oc v . 
We point out that this result is remarkably similar to 
the average value of the lo w energy spectral slopes seen 
in large samples of GRBs ( |Kaneko et al. ||2006j |2008[ ). 
We thus may obtain a natural explanation to this obser- 
vational result, in a model that considers emission from 
the photosphere, once the full spatial scattering positions 
and scattering angles are taken into account (of course, 
with several limitations which are discussed below). 

This paper is organized as follows. In [J2] we briefly 
discuss the properties of the photosphere, and the char- 
acteristic time scales up to which thermal emission is 
expected. In Splwe describe the construction of the prob- 
ability densit ytunction. Th ese sections closely follow the 
treatment by Pe'er ( 2008 1 , and are g iven here for com- 
pleteness. We then calculate in §3.1| the observed spec- 
trum at late times, and show that the energy spectrum 
can be approximated as F v {t) oc v°t~ 2 . We compare 
the analytical predictions to the numerical results in <HJ 
We then discuss the implications of our results on the 
observed GRB prompt emission spectra in f|5j 

2. BASIC CONSIDERATIONS: PHOTOSPHERIC RADIUS 
AND CHARACTERISTIC TIME SCALES IN 
RELATIVISTICALLY EXPANDING PLASMA WIND 

Consider the ejection of a spherically symmetric 
plasma wind from a progenitor characterized by con- 
stant mass loss rate M, that expands with time inde- 
pendent velocity v — f3c. The ejection begins at t = 
from radius r = 0, thus at time t the plasma outer edge 
is at radius r out (t) — f3ct from the center. However, 
here we assume that the plasma wind occupies the en- 
tire space, i.e., r out (t) — > oq_| For constant M and V, 
at r < r out the comoving plasma density is given by 
n'(r) = M/(4irm p vTr 2 ), where r = (1 - (3 2 )- 1 / 2 . We 
further assume that emission of photons occurs deep in- 
side the flow where the optical depth is r 3> 1, as a 
result of unspecified radiative processes. The emitted 
photons are coupled to the flow (e.g., via Compton scat- 
tering), and are assumed to thermalize before escaping 
the plasma once the optical depth becomes low enough. 

4 This assum ption has very little effect on the obtained results; 
see discussion in |Pe'er | ( |2008} . 
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Und er these assumptions, it was shown by |Pe'er | the j et opening angle in G RB outflow, < 9j ~ 0.1 rad 
(2008]), that the optical depth of a photon emitted at 
radius r, angle to the line of sight 9 and propagates to- 
ward the observer is 



r(r,0) 
Here, 



Rd 



sin(0) 



Rd 

2nr 



1 

f2 



Rd 



Mcj T 
4m p /3c ' 



(1) 



(2) 



m p is the proton rest mass and a? is Thomson cross 
section. The last equality in equation [T] holds for T 3> 1 
and small angle to the line of sight, 9 -C tt/2, which 
allows the expansion sin(0) ~ 9 — 9 3 /6. 

The photospheric radius is obtained by setting 
r(r ph> 0) = 1, 



r ph (9) 



Rd 

2^ 



1 

fa 



(3) 



Equation[3]implies that for small viewing angle, 9 <C T 1 
the photospheric radius is angle independent, r p h ~ 
i?d/27rr 2 , while for large angles 9 ^> the photo- 

spheric radius is r p h(9) ~ Rd9 2 /6n. 

Characteristic observed times. Below the photosphere, 
the photons are coupled to the flow, therefore their effec- 
tive propagation velocity in the radial direction is similar 
to the outflow velocity, w /3c. Assuming that a photon is 
emitted at t = 0, r = 0, it decouples the plasma at time 
t = 7"//3cjj Consider photons that propagate towards the 
observer at angle to the line of sight 9. These photons 
are observed at a time delay with respect to a hypothet- 
ical photon that was emitted at t = 0, r = and did not 
suffer any time delay ("trigger" photon), which is given 
by 

At° b -(r,9) = ^-x [l-/3cos(0)]. (4) 

For relativistic outflows, r> 1, thermal photons emit- 
ted from the photospheric radius r p hi&) on the line of 
sight (9 = 0), are seen at a very short time delay with 
respect to the trigger photon, which is given by 

At° b - (r ph , = 0) = ~ ~ 10- 2 LI7 5 s. (5) 

Here, M = L/Tc 2 , and typical parameters characterizing 
GRBs, L = lO^Lergs" 1 and T = 100r 2 were used. 

On the other hand, due to the strong angular depen- 
dence of the photospheric radius, thermal photons emit- 
ted from the photosphere with high angles to the line 
of sight, S > r 1 (and 9 -C 1) are observed at a much 
longer time delay, 



t ob.( 6 >r -l) 



R d (9 



37T/?C 



ZOLT^ei-, s, 



(6) 



where 9 — O.10_i rad. The time scale derived on the 
right hand side of equation [6] is based on the estimate of 

5 This result heavily relies on the assumption of constant outflow 
velocity below the photosphere. In GRBs as well as other astro- 
nomical objects, acceleration episode is expected, which changes 
the characteristic time scale of photon emergence. Nonetheless, we 
use this simplified assumption here, and further discuss it in j|5| 



(e.g., Berger et al. 2003). Equation [6j therefore shows 
that in a relativistically expanding wind with parame- 
ters that can characterize GRBs, thermal emission is ex- 
pected up to tens of seconds following the decay of the 
inner engine. 

3. USE OF PROBABILITY DENSITY FUNCTION IN 
CALCULATING THE EMISSION PROPERTIES 

Any attempt of describing the photospheric emission 
must consider the fact that photons have a finite proba- 
bility of being scattered at every point in space in which 
electrons exist (i.e., in the entire sp ace, and not o nly on 



the photospheric surface). This led Pe'er (2008) to in- 



troduce the concept of a probability density function, as 
a mathematical tool which is needed in calculating the 
observed properties of the photospheric emission, for ex- 
ample the temporal decay of the temperature and flux. 
Since, as we show here, this is a fundamental concept 
which is necessary in an analysis of the properties of the 
photospheric emission, and in particular calculation of 
the expected spectrum, we briefly repeat in this section 
the basic definition and u se o f this function, before cal- 
culating the spectrum in |3.1| 

The thermal photons are advected with the flow below 
the photosphere, until the last scattering event (the de- 
coupling) takes place, at some radius r. For every radius 
r there is an associated probability that the last scat- 
tering event occurs at that particular radius (see below). 
During the last scattering event, a photon is scattered 
into angle 9. An underlying assumption is that photons 
that decouple from the plasma at radius r and scattered 
into angle 9 arc observed at a delay given by equation 
|4j This assumption implies that: (a) the delay time of 
a photon is solely determined by two parameters, the 
last scattering radius and scattering angle (for constant 
outflow velocity); and (b) at any given instance, an ob- 
server sees simultaneously photons emitted from a range 
of radii and angles, all fulfilling the requirement set by 
equation [4] 

The ^-dependence of the photospheric radius implies 
that the probability of a photon to be scattered into 
angle 9 depends on the radius at which the scatter- 
ing event takes place (or vice versa). However, here 
we assume that the probabilities are independent, i.e., 
P(r, 9) = P(r) x P(9). This assumption is made in order 
to simplify the calculation, and is tested against the nu- 
merical results (see ^4] below). While clearly this approx- 
imation has only a limited validity, the results obtained 
are in good agreement with the precise calculation done 
numerically. We further discuss this approximation, as 
well as its limitations in <|5] below. 

Equation [T] implies that the optical depth to scatter- 
ing depends on the radius as r(r) cx r _1 . This optical 
depth is the integral over the probability of a photon 
propagating from radius r to +oo to be scattered, i.e., 
r(r) = J r [dr / dr)dr , from which it is readily found that 
(dr/dr)\ r oc r~ 2 . As a photon propagates in the radial 
direction from radius r to r + Sr, the optical depth in 
the plasma changes by St = (dr / dr)\ r 5r . Therefore, the 
probability of a photon to be scattered as it propagates 
from radius r to r + Sr is given by 



Psc.( r -- r + <5 r ) = 1 — e 



St oc 



Sr 



(7) 
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For the last scattering event to take place at r..r + Sr. 
it is required that the photon does not undergo any ad- 
ditional scattering before it reaches the observer. The 
probability that no additional scattering occurs from ra- 
dius r to the observer is given by exp(— r[r]). The prob- 
ability density function P(r) for the last scattering event 
to occur at radius r, is therefore written as 



P(r) 
P(r) 



= r je~ (ro/r) . 



(8) 



The function P(r) in equation [8] is normalized, 
J °° P(r)dr = 1. Comparison to equation [l] gives the 
proportionality constant, r$ = r p h{9 = 0) = Rd/2irT 2 . 

The probability of a photon to be scattered into angle 9 
is calculated assuming isotropic scattering in the comov- 
ing frame, i.e., da/dil' = Constrj The comoving spatial 
angle is dfl' = sin O'dO'dcf/, and therefore the probability 
of a photon to be scattered to angle 9' (in the comoving 
frame) is dP/d9' oc sin#'. The proportionality constant 
is obtained by integrating over the range < 9' < ir, and 
is equal to 1/2. Thus, the isotropic scattering approxi- 
mation leads to P{9') = (sin 6')/2. 

Assuming that on the average, photons propagate in 
the radial direction, by making Lorentz transformation 
to the observer frame one obtains the probability of scat- 
tering into angle 9 with respect to the flow direction. 
This angle is equal to the observed angle to the line of 
sight, and is given by 

p(6) = p(e') — dcosd ' dcos9 = 

1 ; V ; dcos6" dcos9 d9 2T 2 (l 
Defining u=l — f3 cos 6, equation [9] becomes 

P ^=2TW (10) 
< u < 1 + /?, and the function P(u) in 



sm( 



>cos( 



Note that 1 
equation 



10 



is normalized, P(u)du 



3.1. Spectrum and decay law of the thermal flux at late 

times 

As long as the radiative processes that produce the 
thermal photons deep inside the flow are active, the ob- 
served thermal radiation is dominated by photons emit- 
ted on the line of sight towards the observer. Once 
these radiative processes are terminated, the radiation 
becomes dominated by photons emitted off axis and from 
larger radii, which determine the late time behavior of 
the spectrum and flux. Therefore, the limiting case of 
a (5-function injection, both in time and radius (t = 0, 
r = 0) is expected to closely describe the late time be- 
haviour of the thermal spectrum. We calculate here the 
observed spectrum and flux of the thermal emission at 
late times, under these assumptions. 

Denote by T'(r) the photon comoving temperature, 
its observed temperature is T ob ~ T'(r)T>, where T> = 
[r(l - /3COS0)]- 1 = (ru)- 1 is the Doppler factor. The 
observed photon temperature therefore depends on the 
viewing angle as well as on the radius of decoupling. Be- 
low the photospheric radius, photons lose their energy 

6 This assumption neglects the dipole approximation, and is 
checked numerically to be valid. 



adiaba t ically, hence T'(r) 



oc r 



-2/3 



It was shown by 



Pe'er (2008), that a similar decay law for the photon 



temperature exists even if the energy density in the pho- 
ton field is much smaller than the energy density in the 
electrons (which can in principle be non-relativistic in the 
comoving frame, hence have a different temperature de- 
cay law), resulting from the spatial 3-d expansion of the 
plasma. The adiabatic losses take place only as long as 
the photons propagate at radii smaller than ~ few x rrj. 
While photons that propagate at high angles decouple 
the plasma at much larger radii than r (see eq. pi), 
above few x r , the number of scattering is small, and 
hence the photons maintain their energy. As here we 
are interested in the late time evolution of the spectrum 
and flux, where late time imply t 3> tiy, we can safely 
assume that the comoving temperature of photons that 
dominate the flux at late times is (on the average) con- 
stant, T'(r) = Const. 

Assume that Ao photons are emitted instantaneously 
(a ^-function injection in time) at the center of the ex- 
panding plasma. Each photon is advected with the flow, 
until the last scattering event takes place at radius r and 
into angle 9, after which it propagates freely. The ob- 
served flux density (or differential energy flux) from a 
source at luminosity distance d^ is therefore given by 



ob\ 



_ d 2 F _ hN 
dtdv 



J P(r)dr J P(u)duT° b (r,u) 



T'(r) 

r« 



(11) 

Here, F is the total fluence, and h is Planck's constant; 
the observed frequency v corresponds to the observed 
temperature, v — T ob /h. 

As discussed above, at late times, t ^> ijv, one can 
write T ob = T^/Tu, where Tq is r-independent. Us- 
ing P(r) and P(u) from equations [8] and 10 and the 
identities 5{t ob = ru/(3c) = 5(r = Pct ob /u) x /3c/u and 



S(T° b = T^/Tu) 
[Til becomes 

F„(t° b ^>t N ) 



= 5(u = Tq/TT ) x Tu 2 m, equation 



x 



hN r rg 

hN c f ro 
T 2 J r 2 ' 

_ fict 



Pet 



2/3r 2 « : 



-du 



xS (r 

fejVpc 
87Td'i r 



~ dr 



\r/3ct° b J 



r 2 T° 

T ,2 



H2) 

This equation can further be simplified by using the defi- 
nition of 1 7v from equation [5] and noting that at any given 
instance, the maximum observed temperature is given by 
C = TUTu min = T'/T(l - p) ~ 2r^. Using these in 
equation |T2] leads to the final form, 



F v {t ob » t N ) 



hN c 



2ird 2 L r 



rt N 

t ob 



(13) 



Jmax/' 1 - 



where v max 

Equation [13] is the key finding of this paper. It im- 
plies that atlate times, t ob ^> t^, for a wide frequency 
range v max (tM /t ob ) < v < v maxi the exponent is close 
to unity, and therefore a flat energy spectrum F v cx i/° 
is expected below the thermal peak. This spectrum re- 
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suits from a simultaneous observation of thermal photons 
emitted from a large range of radii and angles to the line 
of sight. We emphasis again that this is a purely geo- 
metrical effect, as no additional radiative proc esses are 
considered. We further note that equation |13| provides, 
in addition, the decay law of the energy flux at late times, 

F v (t ob ^> t N ) <x t ob ~ 2 . 

4. NUMERICAL CALCULATION OF THE FLUX AND 
TEMPERATURE DECAY AT LATE TIMES 

The analytical calculations presented above were 
checked with a numerical code. The code is a Monte- 
Carlo simulation, based on earlier code developed for the 
study of ph oton propagation in relativistically expand - 
ing plasma ( |Pe'er fc Waxman |2004| |Pe'er et al. |2006b| ) . 
This code is essentially identical to the code used fo r 



the numerical calculations that appear in IPe'er I p008 



and a description of it appears there. We give here only 
a basic description of the code, for completeness, before 
presenting the numerical results and a comparison to the 
analytical approximation presented above. 

The code is essentially a Monte-Carlo simulation of 
Compton scattering between photons and electrons. The 
uniqueness of it lies in the fact that it calculates the inter- 
actions during a relativistic, three-dimension expansion 
of the plasma. As a result, the probability of a photon 
to be scattered at any given instance, which is translated 
to the distance traveled by a photon between two con- 
secutive scattering events, depends on the instantaneous 
radius and propagation direction of the photon. In every 
interaction, the full Klein-Nishina cross section is used in 
calculating the outgoing photon energy and propagation 
direction. Since a scattering event is calculated in the 
electron's rest frame, before and after every scattering, 
the photon 4-vector is being Lorentz transformed twice. 
First, into the (local) bulk motion rest frame of the flow 
which assumes an expansion at constant Lorentz factor T 
in the radial direction. A second Lorentz transformation 
is made into the electrons rest frame: in the bulk frame, 
the electrons assume a random velocity direction, with 
velocity drawn from a Maxwellian distribution with tem- 
perature T' el (r) oc r -2 / 3 . The proportionality constant is 
determined using parameters characterizing the prompt 
emission in GRB's. 

Initially, photons are injected into the plasma in a 
random position on the surface of a sphere at radius 
Tinj = i?rf/(2-7rr 2 d). The depth d is taken as d = 20 
in order to ensure that the probability of a photon to es- 
cape without being scattered is smaller than exp(— 20), 
i.e., negligibleP] The initial photon propagation direction 
is random, and its (local) comoving energy at the injec- 
tion radius is equal to the plasma comoving temperature 
at that radius. 

4.1. Numerical results 

We show in figure [l] the positions of the last scattering 
events for N = 3 x 10 6 photons in the r — 9 plane. The 
three contour lines (thin black) are added to the plot in 
order to help demonstrating the probabilities of photons 
to be scattered at a given radius and angle. The thick 



(green) line, is the photospheric radius, as calculated in 
equation [3j Clearly, the photospheric radius provides 
only a first order approximation to the last scattering 
events positions. It is obvious from the figure that pho- 
tons decouple from the plasma at a range of radii and 
angles, necessitating the use of the probability density 
functions introduced in Sj3] We further added to the plot 
(dashed, blue lines), three equal arrival times contours 
(see eq. , for three different values of the observed time. 
These contours demonstrate one of the key results: while 
at early times, t ob < few x tjy, the emission is dominated 
by photons emitted at angles smaller than L _1 , at later 
times, most of the contribution to the emission is from 
photons emitted from higher angles, 6 > F^ 1 . This is 
the origin of the emerging power-law spectrum. 

In preparing the plot, we used parameters that can 
characterized GRBs, such as initial expansion radius — 
10 8 cm, luminosity L — 10 52 ergs _1 and Lorentz factor 
L = 400. In accordance to equation [3j the photospheric 
radius is angle-independent at angles smaller than = 
2.5 x 10 -3 rad, and is oc 2 at higher angles. We point 
out that in figure [T] the radius of the last scattering event 
is normalized to r — R^/2ttT 2 (r = 9 x 10 10 cm for 
the parameters chosen here). Therefore, results obtained 
for arbitrary values of the free model parameters (L and 
F) that characterize astrophysical transients other than 
GRBs, such as AGNs or microquasars are similar to the 
ones presented. 

In figure [2] we present the resulting observed spectrum, 
F v calculated at different times. The solid lines show 
the full numerical results at observed times t ob /tjv = 
10, 10 2 ,10 3 ,10 4 . The dashed (thin) lines are the ana- 
lytical approximation given by equation |13| In calcu- 
lating the maximum observed temperature, T° b ax , we 
used the standard "fireball" model to determine the 
temperature evolution below t he photosphere (see, e.g., 

Thus, the plasma as- 
and the saturation ra- 



7 In fact, it was shown by |Pe'er |f2008| ) that the average number 
of scattering prior to photon escape is « 2d. 



Piran ||2005[ |MeszarosJ |2006j) . 
sumes to accelerate between 
dius r s = Vri, and continues with constant outflow ve- 
locity at larger radii. For the values of the parameters 
chosen, the plasma comoving temperature at the photo- 
spheric radius r is T'(r ) = (i/47rr 2 ca) 1 / 4 (r /r s ) _2 / 3 = 
0.54L" 5/12 r l 1/ 8 6 (r/400) 5 / 3 keV, and the maximum ex- 
pected temperature is T£ h ax = 2rT'(r ) = 430 keV. 

The numerical results are indeed in very good agree- 
ment with the analytical approximation. We note that 
for the parameters chosen here, which can characterize 
many GRBs, % s» 10~ 5 s, and therefore the observed 
characteristic times of < second (in a single pulse) are 
translated into t ob /ijy < 10 4 . At early times, the ex- 
ponential decay occurs at (relatively) high energies, and 
thus a flat spectrum is not seen. However, at later times, 
the spectrum becomes flat over a wide energy range, con- 
sistent with the analytical approximations. 

At intermediate times there is a slight discrepancy be- 
tween the numerical results and the analytical approxi- 
mation, both at the high and low energies. This discrep- 
ancy demonstrates the li mita tion of the delta-function 
approximation used in eq. 13] At high energies, the nu- 
merical results show "bending" at early times, which is 
not captured by the analytical approximation. This re- 
sults from the the (5-function approximation in time used 
(eq. [5]). The underlying assumption used in deriving 
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Angle to the line of sight, > 

Fig. 1. — Position of the last scattering event in 8-r plane for 
10 6 ' 5 events. Parameters chosen in producing the plot are L = 
10 52 erg s — 1 and T = 400. The solid (green) line is the photospheric 
radius, calculated in equation[3] Clearly, the last scattering events 
occur over a wide range of radii and angles. The photospheric 
radius gives only a first order approximation to the position of 
these events. The contour lines are added to the plot in order to 
indicate the density of the emerging photons radii and angles. The 
blue (dashed) lines are equal arrival time contours, calculated in 
equation[5] While at early times, t > tjv the emission is dominated 
by photons scattered at angles 9 < r — 1 , at late times, photons 
emitted from large angles dominate the flux, a fact that gives rise 
to the flattening of the spectrum. 

equation [5] is that the observed delay time is solely de- 
termined by the radius and angle of the last scattering 
event. In reality, of course, as the photons diffuse below 
the photosphere, the time delay is determined by the full 
history of the photon propagation below the photosphere. 
At high energies, the spectrum is dominated by photons 
emitted with high Dopplcr shift, i.e., at small viewing 
angles; as at late times there are relatively very few such 
photons, dispersion in the last scattering position and 
angle leads to a high energy decay. 

At low energies, the discrepancy between the analytical 
prediction and the numerical results is explained due to a 
combination of two phenomena. First, the analytical ap- 
proximation used does not consider the coupling between 
the probability density functions P{r) and P(9). As is 
clear from figure [l] for any given radius r, the observed 
viewing angle 9 is limited from above. A second source of 
spread lies in the fact that the comoving photons energy 
distribution also have an internal spread, (although, the 
average comoving energy is constant above few x ro). 
As a result of these two effects, the low energy spectrum 
is also slightly bended with respect to the simple power 
law prediction of the analytical approximation. For a full 
treatment of these effects, one needs to obtain a full so- 
lution of the diffusion equation, which will be presented 
elsewhere. 

In spite of being only a first order approximation, it is 
clear from figure [2] that the analytical predictions are 
in very good agreement with the numerical results of 
the observed spectrum. At intermediate times, the spec- 
trum approaches a power law distribution F v cx is over 
a relatively wide energy band. Even at very early times, 
t > ijv, it is clear that the observed spectrum deviates 
significantly from the classical Planck function. 

5. SUMMARY AND DISCUSSION 

In this manuscript, we have studied both analytically 
and numerically the spectrum resulting from a photo- 




10 10 
E ob =h v [eV] 

Fig. 2. — Observed spectrum at different times, t/tfj = 
10, 10 2 , fO 3 , 10 4 . The solid lines are the results of the numeri- 
cal simulation, and the dashed ( thin ) lines are the analytical ap- 
proximation derived in equation |13| For initial expansion radius 
rj = 10 8 cm and other parameters same as in figure [l] as consid- 
ered here, T^j^ = 430 keV (see text for further details). At late 
times, ( > if;, the spectrum becomes flat, F v oc u° over a wide 
spectral range below T^ x . This spectrum is very different than 
the Rayleigh-Jeans tail of a thermal spectrum. 

spheric emission in relativistically expanding plasma. We 
showed that even for purely thermal distribution of the 
comoving photon spectrum, aberration of the light re- 
sults in an observed spectrum at late times that is very 
different from a Planck distribution. We showed that at 
late times, t 3> tfj, the observed spectrum approaches 
a power law, with power law index F v oc v over a 
wide spectral range. This result is remarkably similar 
to the observed low energy spectral index of ma ny GRBs 
( |Preece et a/~|[2000j |Kaneko et aT]^00E\ [2008] ) . 

The origin of this non-trivial result lies m two facts: 
First, at late times, the emission is dominated by pho- 
tons emitted off- axis (from angles 9 ^> r _1 , see figure 
[I]). Due to the strong dependence of the Doppler shift 
on the angle to the line of sight, these photons are seen 
at much lower energies than photons emitted on-axis. 
Second, at any given instance, an observer sees simulta- 
neously photons that are emitted from a range of radii 
and angles (see eq. J5| . As each photon has a finite prob- 
ability of being emitted from a given radius, and is seen 
at a particular Doppler shift, the observed spectrum and 
flux can only be described in terms of probability den- 
sity functions (see These functions provide a math- 
ematical tool to describe the probability of photons to 
be emitted from radius r and into angle 9. Using these 
functions, we calculated an an alytical approximation to 
the observed spectrum (eq. 13), which is the main result 
of this work. The analytical approximation was tested 
with a Monte-Carlo simulation that tracks the evolution 
of thermal photons in relativistically expanding plasma 
(Q. The simplified analytical calculations are found to 
be in very good agreement with the accurate numerical 
results (see figure [5]). 

In spite of the success in reproducing the low energy 
spectral slopes at late times, the theory is still not com- 
pleted. For parameters characterizing GRBs, the ther- 
mal peak ( T^ax) naturally falls at the sub-MeV range 
(see [4.1), hence the observed peak can naturally be ex- 
plained as having a thermal origin. However, a noticeable 
drawback of the theory as stated in this manuscript, is 
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that the same parameters lead to very short character- 
istic time scale, ijv ~ 10~ 5 s. Since jtbove t^ the flux 
decays rapidly, F u (t) cx t~ 2 



(see eq. 



13 ), one expects a 
relatively weak thermal signal at f > t^, which can still 
be very short time. We note though, that the calcula- 
tion of ij\r in equation [5] is based on the assumption of 
constant outflow velocity. This assumption is too simpli- 
fied: in order to reach high Lorentz factor, the plasma 
needs to undergo an acceleration phase, and so the av- 
erage Lorentz factor below the photosphere is less than 
the terminal Lorentz factor. As a result, we expect that 
in practice the characteristic time scale relevant for ther- 
mal emission in GRBs is longer than the one considered 
in equation [5j 

The exact delay time depends on several uncertain 
conditions. One is the content of the fireball: for ex- 
ample, Poynting-flux dominated fireball is expected to 
have slower acceleration than matter dominated fireball 
jPrenkhahn ]|2002l [Drenkhahn fc Spruit ||2002] jGiannios 



fe Spruit ||2005||2066| ). Another is the baryon load ( |loka 



the jet dMorsony et al. 2007 


Lazzati et al. 12009 


Mizuta 


et al. ||20IUp: although m this manuscript we considered 



terized by regions of higher and lower densities, and is 
thus not steady (edge effects due to the finite opening 
angle may also play a role at late times). 

An additional source of discrepancy between the the- 
oretical predictions developed in this paper and the ob- 
served spectrum, lies in the fact that the theory here 
does not consider any additional, non-thermal radiative 
processes. As shown here, photospheric emission is ca- 
pable of reproducing the peak energy and the low en- 
ergy spectral slope (a in the "Band" function) seen in 
GRBs (at late times). However, photospheric emission is 
not capable of of producing high energy photons (above 
T max < MeV), as are seen in some GRBs by the LAT de- 
tector on board the Fermi satellite. The inclusion of high 
energy, non-thermal photons, necessitates additional ra- 
diative mechanisms, that must take place following dissi- 
pation processes that occur above the photosphere (e.g., 
in GRB080916C analysis of the high en ergy emission im- 
ply P oynting dominated outflow; see |Zhang fe Pe'er | 
20091. Additional radiative mechanisms, such as syn- 



chrotron emission or Compton scattering, naturally pro 
duce a broad band energy spectrum, and thus may con- 
tribute not only to the high energy spectrum but to the 
low energy part (below the thermal peak) as well. The 
overall observed spectra below the thermal peak is thus 
generally expected to be hybrid- i.e., composed of both 
thermal and non-thermal parts. The exact contribution 
of the non-thermal part may vary from burst to burst, 
depending on the values of the free model parameters, 
such as the radius of the photosphere, the dissipation 
radius, the strength of the magnetic field, etc. 

In principle, the inclusion of thermal photons con- 
tributes to the high energy, non thermal part of the spec- 
trum as well, as these photons serve as seed photons for 
Compton scattering by the non-thermal electrons and 



hot pairs (Recs & Mcszaros 2005 


Pe'er et al. 12005 2006 


Lazzati & Bcgclman ||2010| |Beloborodov |2U10). The ex- 



act contribution of the thermal photons depend on the 
optical depth at the dissipation radius. If the dissipation 



radius is close to the photosphere, the resulting spectrum 
has a complex shape, that is very different than either a 
thermal spectrum or the optically thin synchrotr on - syn 



al. 2005 



2006| ). On the other hand, if the dissipation 
occurs at large radii, the two components, thermal and 
non-thermal, can be decoupled, a fact that can lead to 
a clear identification of t he thermal component, as is in 
the c ase of GRB090902B ( |Ryde et al. pOiol |Pe'er et al. ' 



chrotron s elf Compton (SSC) model predictions (Pe'er et 

T 



2010) 



One clear prediction of the results presented here, is 
that if a thermal component contributed significantly 
to the observed spectrum, than the low energy spectral 
slope (a) is expected to vary with time: at early times, 
when the inner engine is active, a is expected to be steep 
(close to the Rayleigh- Jeans tail), while at later times, 
F v oc v°. This result is expected to be correlated with a 



possibly observed in several bursts (e.g., Crider et al. 


1997; Ghirlanda et al. |2007| 


Ryde & Pe'er 12009||Page et 


al. |2009p. We note though t 


lat in many GRBs, the data 



the outflows in GRBs are generally not smooth, but very 
fluctuative. Therefore, photons originating from differ- 
ent emission episodes (each characterized by its own mass 
loss ejection rate and its own Lorentz factor) are often 
superimposed. Thus, in fact, thermal photons may only 
be identified using detailed time resolved spectral analy- 
sis of separa te flares. This task was carried by |Ryde te| 



Pe'er ( |2009[ ) , which were indeed able to identify the ther- 
mal component in a large sample of bursts. Once done, 
this analysis method can be fu rther used to dedu ce the 
properties of the GRB outflow ( |Pe'er et q?~||2007| ). 



In addition to the prompt emission phase m GRBs, 
thermal activity may occur as part of the flaring activ- 
ity observed in the ear l y afterglow phase of many GRBs 
(iBurrows et al. ||2005l (Falcone et al. ||2007[ |Chincarini 



et al. ||2010| |Marguttf et al. ||2010| ). 



he exact nature 

ot these hares is currently not yet clear. As it is plausi- 
ble that the flares result from renewed emission from the 
inner core, a renewed thermal emission may occur. The 
analysis presented here may therefore apply in the study 
of the late time flares as well. 

While analysis of GRB prompt emission spectrum 
serves as a major motivation to this work, we note that 
the analysis presented here is general. With modified 
parameters values, the analysis apply to emission from 
any astronomical transient, which is characterized by rel- 
ativistic outflow speeds and high density core. Thus, we 
expect the analysis carried here to be valid for objects 
such as AGNs and microquasars. This is valid because 
the exact nature of the radiative process that produces 
the thermal photons is of no importance, as long as the 
emission occurs deep enough in the flow, where the opti- 
cal depth is high, t ^> 1. 
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